Loading packages

library(ggplot2)
library(dplyr)
library(tidyr)
library(magrittr)
library(purrr)
library(tibble)
library(readr)
library(here)
library(ggrastr)
library(cowplot)
library(princurve)
library(scico)
library(ggridges)

# parameter
print("parameter input:")
## [1] "parameter input:"
print(params$data)
## [1] "data/processed/morphology/umap_absolute_all_drugs_sampled.Rds"
print(params$sample)
## [1] "data/processed/morphology/umap_absolute_all_drugs_tidy_Paclitaxel.Rds"

loading input data and annotation. Note that on the central cluster, with access to the complete data table, the definition of the input can easily be changed. For remote work, the subsampled dataset “umap_drugs_sampled.Rds” is the default choice.

# I wish I could solve my path problems with the here() package, but experienced unreliable behavior 
# PATH = "/dkfz/groups/shared/OE0049/B110-Isilon2/promise/"
PATH = paste0(here::here(), "/")
here::here()
## [1] "/home/rstudio/promise"
#umap_df <- read_rds(paste0(PATH, "data/processed/PhenotypeSpectrum/umap_absolute_all_drugs_tidy.Rds"))
#umap_df <- read_rds(paste0(PATH, "data/processed/PhenotypeSpectrum/umap_absolute_all_drugs_sampled.Rds"))
umap_df <- read_rds(here::here(params$data))
umap_df_sample <- read_rds(here::here(params$sample))

organoid_morphology <- read_delim(here::here("references/imaging/visual_classification_organoids.csv"), ";", escape_double = FALSE, trim_ws = TRUE) %>% 
  dplyr::select(line = organoid, morphology = visual_inspection_v2)

# reading PCA variance information
#organoid_viability <- read_rds(here::here("data/processed/ldc_viability.rds"))
pca_variance <- readRDS(here::here("data/processed/morphology/pca_variance.Rds"))

Data cube

umap_df %>% dplyr::distinct(drug, line) %>% dplyr::count(line)
## # A tibble: 12 x 2
##    line        n
##    <chr>   <int>
##  1 D004T01   528
##  2 D007T01   528
##  3 D010T01   528
##  4 D013T01   528
##  5 D018T01   528
##  6 D019T01   528
##  7 D020T01   528
##  8 D020T02   528
##  9 D022T01   528
## 10 D027T01   528
## 11 D030T01   528
## 12 D046T01   528

PCA

pca_variance %>% 
  head(50) %>%
  ggplot(aes(PC, sum_explained)) + 
  geom_point() + 
  geom_hline(yintercept = 0.8, linetype = "dashed") + 
  geom_vline(xintercept = 25) +
  theme_classic() + 
  ggsave(here("reports/figures", params$output, "pca_var.pdf"), width = 4, height = 4)

Partition inspection

We are able to observe 4 partitions in our data. After manual inspection, it becomes cleat that the two smallest partitions are mostly consisting of

umap_df %>% 
  ggplot(aes(v1, v2, color = factor(partition))) + 
  geom_point_rast(alpha = 0.5, size = 0.35) + 
  scale_color_brewer(type = "qual", palette = "Set2") +
  theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       color = "partition") + 
  theme(legend.position = "bottom") + 
    coord_fixed()

umap_df %>% 
  dplyr::count(partition) %>% 
  mutate(ratio = n/sum(n)) %>% 
  arrange(desc(ratio))
## # A tibble: 1 x 3
##   partition      n ratio
##   <fct>      <int> <dbl>
## 1 1         271308     1

I remove 2 partitions from all main figures for ease of reading. Below, it is easy to toggle the removal of partitions on and off to make sure this filtering step is robust

gg_cluster <- umap_df %>%
  filter(partition %in% c(1,2)) %>%
  ggplot(aes(v1, v2, color = factor(cluster))) + 
  ggrastr::geom_point_rast(alpha = 0.5, size = 0.35) + 
  scale_color_manual(values = c(RColorBrewer::brewer.pal(12, "Set3"), "#fb9a99")) +
  cowplot::theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       color = "partition") + 
  theme(legend.position = "bottom") + 
    coord_fixed()

gg_cluster + ggsave(here("reports/figures", params$output, "gg_cluster.pdf"), width = 4, height = 4)

Organoid Size Distributions

I plot a size-distribution.

gg_size_dist <- umap_df %>% 
  filter(partition %in% c(1,2)) %>%
  ggplot(aes(size)) + 
  geom_histogram() + 
  theme_cowplot() + 
  labs(caption = "all treatments, downsampled")

gg_size_dist_log <- gg_size_dist + 
  scale_x_log10()

gg_size_dist_log + 
  ggsave(here("reports/figures", params$output, "gg_size_dist.pdf"), width = 4, height = 4)

I add the eCDF.

df <- umap_df %>% filter(partition %in% c(1,2))

gg_ecdf <- ggplot(df %>% filter(drug == "DMSO")) +
  stat_ecdf(aes(x = size, group = line), 
              geom = "step", size = 1) +
  #scale_color_manual(values = c("#00AFBB", "#E7B800"))+
  labs(y = "f(size)",
       x = "organoid size [pixels]") + 
  theme_cowplot()

gg_ecdf

# variance_explained %>% ggplot(aes(PC, var_explained)) + geom_point() + geom_vline(xintercept = 25) + scale_y_sqrt()


# element wise multiplication
mat = components %>% dplyr::select(-X1) %>% as.matrix() %>% .[,1:25] %>% matrix(ncol = 25)
vector = sqrt(variance_explained$eigenvalue) %>% .[1:25]
df = mat * vector

index = df %>% abs() %>% rowSums() %>% tibble(max = ., index = 1:length(.)) %>% arrange(desc(max)) %>% head(25) %>% .$index

pheatmap::pheatmap(mat = annotated_PC %>% dplyr::select(name, PC1:PC25) %>% .[index,] %>% remove_rownames() %>%  tibble::column_to_rownames('name'),
                   annotation_row = annotated_PC %>% dplyr::select(name, class, channel) %>% .[index,] %>% tibble::column_to_rownames('name'),
                   cluster_cols = FALSE)

For more details about distributions, please refer to *reports/Phenotypespectrum/‘xyz’_dist.pdf*.

line_param <- umap_df %>% filter(partition %in% c(1,2)) %>% 
  #filter(drug == "DMSO") %>%
  nest(-line, -replicate) %>% 
  mutate(fit = map(data, ~ fitdistrplus::fitdist(.x$size, "lnorm")),
         param = map(fit, ~ .x$estimate %>% broom::tidy()))

df <- line_param %>% unnest(param) %>% 
  filter(names == "meanlog") %>% 
  group_by(line) %>% 
  mutate(mean_meanlog = mean(x)) %>% 
  arrange(mean_meanlog) %>% 
  ungroup() %>%
  mutate(line = factor(line, levels = .$line %>% unique()))

organoid_size_fit <- df %>% dplyr::select(line, replicate, names, x, mean_meanlog)
# organoid_size_fit %>% saveRDS(here::here("data/processed/morphology/organoid_size.Rds"))
# organoid_size_fit <- readRDS(here::here("data/processed/morphology/organoid_size.Rds"))

organoid_size_factor <- organoid_size_fit$line %>% levels()
df <- organoid_size_fit

df <- df %>% 
    dplyr::select(line, replicate, x) %>%
    # tidyr::pivot_wider(names_from = replicate, 
    #             values_from = x)
    tidyr::spread(key = replicate, value = x)

r_size = df %>% ungroup() %>% dplyr::select(-line) %>% as.matrix %>% cor() %>% min()

gg_size_replicate <- df %>% 
  ggplot(aes(`1`, `2`)) +
  geom_smooth(method = "lm", se = FALSE, color = "grey") +
  geom_point() + 
  ggrepel::geom_text_repel(data = df %>% filter(line %in% c("D021T01", "D046T01")), aes(label = line)) + 
  theme_cowplot() + 
  labs(x = "Replicate 1, ln(size)",
       y = "Replicate 2, ln(size)",
       caption = paste0("all treatments, downsampled, 2 replicates, r= ", round(r_size, 2))) + 
  coord_fixed(ratio = 1) 
  #geom_abline(slope = 1, color = "grey")
  

gg_size_replicate +
  ggsave(here("reports/figures", params$output, "gg_size_correlation.pdf"), width = 4, height = 4)

organoid_size_factor_09 <- umap_df %>% filter(partition %in% c(1,2)) %>% group_by(line) %>% 
  summarise(x = quantile(size_log, 0.9)) %>% 
  #summarise(x = mean(size_log)) %>% 
  arrange(x) %>% .$line


gg_size_dist_morph_ridge <- umap_df_sample %>% filter(partition %in% c(1,2)) %>% filter(drug == "DMSO") %>% 
  mutate(line = factor(line, levels = organoid_size_factor_09)) %>% 
  ggplot() +
  geom_density_ridges_gradient(aes(y = line, x = size_log, fill = stat(x)), scale = 1) +
  #geom_density(aes(x = size_log, group = replicate, color = morphological_class)) + 
  #facet_wrap(~ line) + 
  scale_fill_viridis_c() +
  labs(caption = "DMSO treated organoids",
       x = "ln(size)",
       fill = "size") + 
  theme(legend.position = "bottom") +
  theme_cowplot() + 
  labs(caption = "DMSO treated, downsampled") +
  coord_fixed(ratio = 1)

gg_size_dist_morph_ridge + ggsave(here("reports/figures", params$output, "gg_size_dist_morph_ridge.pdf"), width = 4, height = 4)

gg_size_dist_morph_violin <- umap_df_sample %>% filter(partition %in% c(1,2)) %>% filter(drug == "DMSO") %>% 
  mutate(line = factor(line, levels = organoid_size_factor_09)) %>% 
  ggplot() +
  geom_violin(aes(y = size_log, x = line, fill = stat(x)), scale = 1) +
  #geom_density(aes(x = size_log, group = replicate, color = morphological_class)) + 
  #facet_wrap(~ line) + 
  scale_fill_viridis_c() +
  labs(caption = "DMSO treated organoids",
       x = "ln(size)",
       fill = "size") + 
  theme(legend.position = "bottom") +
  theme_cowplot() + 
  labs(caption = "DMSO treated, downsampled") +
  coord_fixed(ratio = 1)

gg_size_dist_morph_ridge_flip <- umap_df_sample %>% filter(partition %in% c(1,2)) %>% filter(drug == "DMSO") %>% 
  mutate(line = factor(line, levels = organoid_size_factor_09)) %>% 
  ggplot() +
  geom_density_ridges_gradient(aes(y = line, x = size_log, fill = stat(x)), scale = 1) +
  #geom_density(aes(x = size_log, group = replicate, color = morphological_class)) + 
  #facet_wrap(~ line) + 
  scale_fill_viridis_c() +
  labs(caption = "DMSO treated organoids",
       x = "ln(size)",
       fill = "size") + 
  theme(legend.position = "bottom") +
  theme_cowplot() + 
  labs(caption = "DMSO treated, downsampled") +
  coord_fixed(ratio = 1) + 
  coord_flip()

gg_size_dist_morph_ridge_flip + ggsave(here("reports/figures", params$output, "gg_size_dist_morph_ridge_flip.pdf"), width = 8, height = 4)

umap_size <- function(umap){
  umap %>%
  #filter(Size < 1000) %>%
  ggplot(aes(v1, v2, color = size_log)) + 
  geom_point_rast(alpha = 0.5, size = 0.35) + 
  scale_color_viridis_c() +
  theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       color = "ln(size)") + 
  theme(legend.position = "bottom") + 
    coord_fixed()
}

gg_size <- umap_size(umap_df) + 
  labs(caption = "all treatments, downsampled")

gg_size  + ggsave(here("reports/figures", params$output, "gg_size_all.pdf"), width = 4, height = 4)

In general, DMSO treated organoid lines cover the same latent space than drug treated organoids. This is likely influenced by the large number of untreated organoids in the dataset.

df <- umap_df 

gg_size_supp <- df %>%
  mutate(drug = if_else(drug == "DMSO", "DMSO", "other")) %>%
  ggplot(aes(v1, v2, color = size_log)) + 
  geom_point_rast(alpha = 0.5, size = 0.35) + 
  scale_color_viridis_c() +
  theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2") + 
  theme(legend.position = "bottom") +
  facet_wrap(~ drug) + 
  labs(caption = "downsampled")

gg_size_supp + ggsave(paste0(PATH, "reports/figures/imaging/gg_size_treatment.pdf"))

Paclitaxel case study

drug_size <- umap_df_sample %>% filter(partition %in% c(1,2)) %>% filter(drug == "DMSO" | drug == "Paclitaxel") %>% 
  mutate(concentration = ifelse(drug == "DMSO", 0, concentration)) %>%
  #filter(morphological_class == "disorganized") %>% 
  #filter(morphological_class != "other") %>% 
  mutate(concentration = factor(concentration, levels = c("0", "0.0016", "0.008", "0.04", "0.2", "1.0")))

ggdrug_size <- ggplot(drug_size) +
  geom_density(aes(x = log(size), group = concentration, color = concentration), size = 1.5) + 
  scico::scale_color_scico_d() +
  theme_cowplot() + 
  #scale_x_continuous(limits = c(0, 15000)) +
  theme(legend.position = "bottom") + 
  labs(color = "Paclitaxel Concentration Factor",
       title = "Organoid size distribution",
       x = "ln(size)")

ggdrug_size + ggsave(here::here("reports/figures/imaging/size_treatment.pdf"), width = 5, height = 5)

drug_count <- drug_size %>%
  dplyr::count(concentration, line, replicate, well)

ggdrug_count <- ggplot(drug_count) +
  geom_density(aes(x = n, group = concentration, color = concentration), size = 1.5) + 
  scico::scale_color_scico_d() +
  theme_cowplot() + 
  theme(legend.position = "bottom") + 
  labs(color = "Paclitaxel Concentration Factor",
       title = "Organoid count distribution",
       x = "number of objects per well") + 
  scale_x_log10()

ggdrug_count + ggsave(here::here("reports/figures/imaging/count_treatment.pdf"), width = 5, height = 5)

set.seed(234)
loi = c("D022T01", "D046T01")


df <- umap_df_sample  %>%
  filter(partition %in% c(1,2)) %>%
  filter(drug == "DMSO" | drug == "Paclitaxel") %>% 
  mutate(concentration = ifelse(drug == "DMSO", 0, concentration)) %>% 
  filter(line == loi)

gg_drug <- umap_df_sample %>% filter(partition %in% c(1,2)) %>%
  dplyr::select(-line, -concentration) %>%
  ggplot(aes(v1, v2)) + 
  geom_point_rast(alpha = 1, size = 0.35, color = "#f1f1f1") + 
  geom_point_rast(data = df  %>%
    group_by(concentration) %>% 
    sample_n(1000, replace = TRUE),
  aes(color = concentration),alpha = 1, size = 1.5, shape=16) + 
  #facet_wrap( ~ concentration, ncol = 1) + 
  #scale_color_brewer(type = "seq", palette = "YlOrRd") + 
  #geom_density2d(color = "black") + 
  theme_classic() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       caption = paste0(paste(loi, collapse=" "), ", Paclitaxel"),
       color = "Concentration Factor") + 
  scico::scale_color_scico_d() + 
  facet_wrap(~ line, ncol = 2) +
  theme(legend.position = "bottom") +
  #theme_cowplot(font_size = 8) + 
  theme(legend.position = "bottom") + 
  coord_equal()
 
gg_drug + ggsave(here::here("reports/figures/imaging/paclitaxel_shift.pdf"), width = 10, height = 5)

loi <- c("D022T01", "D046T01")

drug_size_param <- drug_size %>% 
  nest(-concentration, -line) %>% 
  mutate(fit = map(data, ~ fitdistrplus::fitdist(.x$size, "lnorm")),
         param = map(fit, ~ .x$estimate %>% broom::tidy()))

df <- drug_size_param %>% unnest(param) %>% 
  filter(names == "meanlog") %>%
  mutate(concentration = factor(concentration, levels = c("0", "0.0016", "0.008", "0.04", "0.2", "1.0"))) 

df <- drug_size %>% 
  group_by(drug, line, concentration) %>% 
  summarise(x = median(size_log))

gg_size_drug <- df %>%
  filter(line %in% loi) %>%
  #mutate(concentration = as.numeric(as.character(concentration))) %>%
  ggplot(aes(concentration, x)) + 
    geom_point(color = "grey") + 
  geom_line(data = df %>% dplyr::rename(line_h = line) , aes(group = line_h), color = "grey") +
  geom_point(data = df %>% dplyr::rename(line_h = line), color = "grey") +
  geom_point(color = "black") +
  geom_line(aes(group = line), color = "black") +
  labs(y = 'mu ln(size)' ,
       x = "Concentration Factor",
       caption = paste0(paste(loi, collapse=" "), ", Paclitaxel")) +
  facet_wrap(~ line)+
  theme_cowplot()

df %>% write_csv(here::here("reports/tables/source_data_paclitaxel_doseresponse_2c.csv"))

gg_size_drug + ggsave(here::here("reports/figures/imaging/paclitaxel_doseresponse.pdf"), width = 10, height = 5)

## organoid viability classifier

umap_ldc <- function(umap){
  umap %>%
  #filter(Size < 1000) %>%
  ggplot(aes(v1, v2, color = prob_dead)) + 
  geom_point_rast(alpha = 0.5, size = 0.35) + 
  scale_color_scico(palette = 'lajolla') + #lajolla #vikO
  theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       color = "p(dead)") + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_colourbar(barwidth = 0.5, barheight = 10))
}

gg_ldc <- umap_ldc(umap_df) +  coord_equal()

gg_ldc + 
  ggsave(here::here("reports/figures/drug_viability/ldc_umap.pdf"), width = 6, height = 6)

umap_df %>% 
  ggplot(aes(size_log, permeability)) + 
  geom_point_rast(alpha = 0.5, size = 0.35) + 
  geom_smooth() + 
  theme_cowplot()

umap_df %>% 
  ggplot(aes(size_log, prob_live)) + 
  #geom_jitter_rast(alpha = 0.5, size = 0.35) + 
  geom_smooth() + 
  theme_cowplot()

df <- umap_df %>%
  group_by(plate, well) %>%
  summarise(actin = mean(actin), 
            permeability = mean(permeability),
            dapi = mean(dapi),
            size = mean(size_log),
            prob_live = mean(prob_live)) 

df %>%
  ggplot(aes(size, prob_live)) + 
  geom_point_rast(alpha = 0.2) +
  geom_smooth()

df %>%
  ggplot(aes(actin, prob_live)) + 
  geom_point_rast(alpha = 0.2)

df %>%
  ggplot(aes(dapi, prob_live)) + 
  geom_point_rast(alpha = 0.2) + 
  geom_smooth()

df %>%
  ggplot(aes(permeability, prob_live)) + 
  geom_point_rast(alpha = 0.2) + 
  geom_smooth()

CTG data LDC data average size data

Organoid Heterogeneity

I plot 2 organoid lines treated with DMSO control

set.seed(234)


df <- umap_df  %>% filter(partition %in% c(1,2)) %>%
  mutate(cystic = if_else(line == "D013T01" &  well == "D24" & plate == "D013T01P001L02", TRUE, FALSE)) %>%
  mutate(compact = if_else(line == "D046T01" &  well == "D24" & plate == "D046T01P007L02", TRUE, FALSE))

gg_cys_comp <- df %>% 
  sample_frac(0.01) %>%
  ggplot(aes(v1, v2, color = size_log)) + 
  #scale_color_brewer(type = "qual", palette = 2) +
  geom_point_rast(alpha = 0.1, size = 0.35) + 
  geom_point_rast(color = "#F4B400", alpha = 1, size = 0.5, data = df %>% filter(cystic == TRUE)) +
  geom_point_rast(color = "#DB4437", alpha = 1, size = 0.5, data = df %>% filter(compact == TRUE)) + 
  scale_color_viridis_c() +
  labs(color = "size",
       caption = "yellow: D013T01, red: D055T01",
       x = "UMAP 1",
       y = "UMAP 2") +
  theme_cowplot() + 
  coord_fixed()

gg_cys_comp

Organoid line differences

I create a single plot showing the two extreme organoid lines and their distribution within the embedding.

set.seed(123)

loi <- c("D019T01", "D007T01",  "D030T01", "D018T01") #c("D055T01", "D007T01",  "D021T01", "D019T01", "D027T01")
loi <- c("D020T01", "D007T01",  "D030T01", "D018T01")
#loi <- umap_df$line %>% unique()

df <- umap_df %>%
  filter(drug == "DMSO") %>% 
  filter(partition %in% c(1,2))


gg_line <- df %>% dplyr::select(-line) %>%
  ggplot(aes(v1, v2)) + 
  geom_point_rast(alpha = 1, size = 0.35, color = "#f1f1f1") + 
  geom_point_rast(data = umap_df %>%
                    filter(drug == "DMSO") %>% 
                   # filter(line %in% c("D021T01")) %>%
    filter(line %in% loi) %>% 
    mutate(line = factor(line, levels = loi)), #%>% 
      #sample_frac(0.1),
    #mutate(line = factor(line, levels = c("D021T01"))),
  aes(color = line),alpha = .4, size = 0.35, shape=16) + 
  facet_wrap( ~ line, ncol =2) +
  scale_color_brewer(type = "qual", palette = "Set2") +
  #scale_color_manual(values = c(c("#D80D12", "#461C01", "#9a4c91", "#70BE6F", "#24345E"))) +   
  #geom_density2d(color = "black") + 
  theme_classic() +
  labs(x = "UMAP 1",
       y = "UMAP 2")+
       #caption = "control treated organoids") + 
  theme_cowplot(font_size = 8) + 
  theme(legend.position = "nothing")

gg_line + ggsave(paste0(PATH, "reports/figures/imaging/gg_size_four_lines.pdf"), width = 4, height = 4)

df <- umap_df %>%
  left_join(organoid_morphology %>% mutate(line = paste0(line, "01"))) %>%
  filter(drug == "DMSO") %>% 
  filter(partition %in% c(1,2))

gg_morph <- df %>% dplyr::select(-line) %>%
  ggplot(aes(v1, v2)) + 
  geom_point_rast(alpha = 1, size = 0.35, color = "#f1f1f1") + 
  geom_point_rast(data = umap_df %>% left_join(organoid_morphology %>% mutate(line = paste0(line, "01"))) %>%
                    filter(drug == "DMSO"), #%>% 
      #sample_frac(0.1),
    #mutate(line = factor(line, levels = c("D021T01"))),
  aes(color = morphology),alpha = .4, size = 0.35, shape=16) + 
  #facet_wrap( ~ line, ncol =2) +
  scale_color_brewer(type = "qual", palette = "Set1") +
  #scale_color_manual(values = c(c("#D80D12", "#461C01", "#9a4c91", "#70BE6F", "#24345E"))) +   
  #geom_density2d(color = "black") + 
  theme_classic() +
  labs(x = "UMAP 1",
       y = "UMAP 2")+
       #caption = "control treated organoids") + 
  theme_cowplot(font_size = 8)

gg_morph + ggsave(paste0(PATH, "reports/figures/imaging/gg_morphology.pdf"), width = 4, height = 4)

df <- umap_df %>%
  left_join(organoid_morphology %>% mutate(line = paste0(line, "01"))) %>%
  filter(drug == "DMSO") %>% 
  filter(partition %in% c(1,2))

gg_screen_id <- df %>% dplyr::select(-line) %>%
  ggplot(aes(v1, v2)) + 
  geom_point_rast(alpha = 1, size = 0.35, color = "#f1f1f1") + 
  geom_point_rast(data = umap_df %>% left_join(organoid_morphology %>% mutate(line = paste0(line, "01"))) %>%
                    filter(drug == "DMSO"), #%>% 
      #sample_frac(0.1),
    #mutate(line = factor(line, levels = c("D021T01"))),
  aes(color = screen_id),alpha = .4, size = 0.35, shape=16) + 
  facet_wrap( ~ line, ncol =4) +
  scale_color_brewer(type = "qual", palette = "Set2") +
  #scale_color_manual(values = c(c("#D80D12", "#461C01", "#9a4c91", "#70BE6F", "#24345E"))) +   
  #geom_density2d(color = "black") + 
  theme_classic() +
  labs(x = "UMAP 1",
       y = "UMAP 2")+
       #caption = "control treated organoids") + 
  theme_cowplot(font_size = 8) + 
  theme(legend.position = "bottom")

gg_screen_id + ggsave(paste0(PATH, "reports/figures/imaging/gg_screen_id.pdf"), width = 4, height = 4)

Dye Intensity

umap_df %>%
  ggplot(aes(actin, fill = screen_id)) +
  geom_density(alpha = 0.25) + 
  scale_fill_brewer(type = 'qual') +
  facet_wrap(~ line) +
  theme_cowplot()

umap_df %>%
  ggplot(aes(permeability, fill = screen_id)) +
  geom_density(alpha = 0.25) + 
  scale_fill_brewer(type = 'qual') +
  facet_wrap(~ line) +
  theme_cowplot()

umap_df %>%
  ggplot(aes(dapi, fill = screen_id)) +
  geom_density(alpha = 0.25) + 
  scale_fill_brewer(type = 'qual') +
  facet_wrap(~ line) +
  theme_cowplot()

set.seed(123)

umap_df %>%
  filter(drug == "DMSO") %>%
  #sample_n(10000) %>%
  ggplot(aes(v1, v2, color = permeability)) + 
  geom_point_rast(alpha = 0.75, size = 0.35) +
   scale_colour_gradient(low = "white", high = "green") + 
  theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       title = "Permeability staining intensity") + 
  theme(legend.position = "bottom") + 
  ggsave(paste0(PATH, "reports/figures/imaging/gg_permeability.pdf"), width = 4, height = 4)

umap_df %>%
  filter(drug == "DMSO") %>%
  #sample_n(10000) %>%
  ggplot(aes(v1, v2, color = actin)) + 
  geom_point_rast(alpha = 0.75, size = 0.35) +
   scale_colour_gradient(low = "white", high = "red") + 
  theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       title = "Actin staining intensity") + 
  theme(legend.position = "bottom") + 
  ggsave(paste0(PATH, "reports/figures/imaging/gg_actin.pdf"), width = 4, height = 4)

umap_df %>%
  filter(drug == "DMSO") %>%
  #sample_n(10000) %>%
  ggplot(aes(v1, v2, color = dapi)) + 
  geom_point_rast(alpha = 0.75, size = 0.35) +
   scale_colour_gradient(low = "white", high = "blue") + 
  theme_cowplot() +
  labs(x = "UMAP 1",
       y = "UMAP 2",
       title = "DNA staining intensity") + 
  theme(legend.position = "bottom") + 
  ggsave(paste0(PATH, "reports/figures/imaging/gg_dapi.pdf"), width = 4, height = 4)

correlation of organoid size, viability and LDC classification

For this analysis I am comparing orgaonid size, CTG data and LDC classification results

utils::data("organoid_viability", package = "SCOPEAnalysis")
auc_ctg <- auc_ctg %>% filter(!Line %in% c("D021T01", "D054T01", "D055T01", "D015T01", "D052T01")) %>% 
  mutate(id = paste0(Line, "_", library, "_", well)) %>%
  arrange((id)) %>%
  group_by(id) %>% 
  summarise(viability = mean(viability)) %>%
  dplyr::select(id, ctg = viability)
mortality <- mortality %>% filter(!Line %in% c("D021T01", "D054T01", "D055T01", "D015T01", "D052T01")) %>% 
  mutate(id = paste0(Line, "_", Layout, "_", Well.ID)) %>%
  group_by(id) %>% 
  summarise(Percent.Live = mean(Percent.Live)) %>%
  dplyr::select(ldc = Percent.Live, id)
organoid_size_drug <- readRDS(here::here("data/processed/morphology/organoid_size_drug.Rds")) %>%
  mutate(id = paste0(line, "_", substr(plate, 12,14), "_", well)) %>% 
  group_by(id) %>% 
  summarise(mean = mean(mean)) %>%
  dplyr::select(size = mean, id)
organoid_intensity_drug <- readRDS(here::here("data/processed/morphology/organoid_intensity_drug.Rds")) %>%
  mutate(id = paste0(line, "_", substr(plate, 12,14), "_", well)) %>% 
  group_by(id) %>% 
  summarise(permeability = mean(permeability),
            dapi = mean(dapi),
            actin = mean(actin)) %>%
  dplyr::select(dapi, actin, permeability, id)

df <- left_join(auc_ctg, mortality) %>% 
  left_join(organoid_size_drug) %>%
  left_join(organoid_intensity_drug)

df %>%
  ggplot(aes(ldc, ctg)) + 
  geom_point() + 
  theme(legend.title = element_blank()) + xlim(c(0, 1)) + ylim(c(0, 1)) + 
  geom_hline(yintercept = 1) + 
  geom_vline(xintercept = 1) +
  theme_cowplot() + 
   coord_equal() + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ctg_ldc_well.pdf"), width = 4, height = 4)

df %>%
  mutate(line = substr(id, 1, 7)) %>%
  ggplot(aes(size, ctg)) + 
  geom_point(color = "grey", alpha = 0.5) + 
  geom_point(aes(color = line), data = df %>% 
               mutate(line = substr(id, 1, 7)) %>% 
               filter(line %in% c("D007T01", "D018T01", "D019T01"))) + 
  geom_smooth(se = FALSE, 
              method = "lm", 
              aes(color = line), data = df %>% 
               mutate(line = substr(id, 1, 7)) %>% 
               filter(line %in% c("D007T01", "D018T01", "D019T01"))) +
  geom_hline(yintercept = 1) + 
  #geom_vline(xintercept = 1) +
  #geom_smooth(se = FALSE) +
  #theme(legend.title = element_blank()) + ylim(c(0, 1)) +
  theme_cowplot() + 
  theme(legend.position = "bottom") +
  coord_fixed(ratio = 2) +
  scale_color_brewer(type = "qual", palette = "Set2") +
  labs(caption = paste0("overall correlation ", cor(df$ctg, df$size, use = "complete.obs") %>% round(digits = 4))) + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ctg_size_well_highlight.pdf"), width = 6, height = 6)

df %>%
  mutate(line = substr(id, 1, 7)) %>%
  ggplot(aes(ldc, ctg)) + 
  geom_point(color = "grey", alpha = 0.5) + 
  geom_point(aes(color = line), data = df %>% 
               mutate(line = substr(id, 1, 7)) %>% 
               filter(line %in% c("D007T01", "D018T01", "D019T01"))) + 
  geom_smooth(se = FALSE, 
              method = "lm", 
              aes(color = line), data = df %>% 
               mutate(line = substr(id, 1, 7)) %>% 
               filter(line %in% c("D007T01", "D018T01", "D019T01"))) +
  geom_hline(yintercept = 1) + 
  geom_vline(xintercept = 1) +
  #geom_smooth(se = FALSE) +
  theme(legend.title = element_blank()) + ylim(c(0, 1)) + xlim(c(0, 1)) +
  theme_cowplot() + 
   coord_equal() + 
  scale_color_brewer(type = "qual", palette = "Set2") +
  labs(caption = paste0("overall correlation ", cor(df$ctg, df$ldc, use = "complete.obs") %>% round(digits = 4))) + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ctg_ldc_well_highlight.pdf"), width = 6, height = 6)

correlation on line-level

df %>%
  mutate(line = substr(id, 1, 7)) %>%
  ggplot(aes(size, ctg)) + 
  geom_point() + 
  geom_smooth(se = FALSE) +
  geom_hline(yintercept = 1) + 
  #geom_vline(xintercept = 1) +
  geom_smooth(se = FALSE) +
  theme(legend.title = element_blank()) + ylim(c(0, 1)) +
  theme_cowplot() + 
  facet_wrap(~line) + 
  labs(caption = paste0("correlation ", cor(df$ctg, df$size, use = "complete.obs") %>% round(digits = 4))) + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ctg_size_well_perline.pdf"), width = 6, height = 6)

df %>%
  mutate(line = substr(id, 1, 7)) %>%
  ggplot(aes(size, ldc)) + 
  geom_point() + 
  geom_hline(yintercept = 1) + 
  geom_smooth(se = FALSE) +
  theme_cowplot() + 
  facet_wrap(~line) + 
  labs(caption = paste0("correlation ", cor(df$ldc, df$size, use = "complete.obs") %>% round(digits = 4))) + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ldc_size_well_perline.pdf"), width = 6, height = 6)

df %>%
  mutate(line = substr(id, 1, 7)) %>%
  ggplot(aes(dapi, ctg)) + 
  geom_point() + 
  geom_hline(yintercept = 1) + 
  geom_smooth(se = FALSE) +
  theme_cowplot() + 
  facet_wrap(~line) + 
  labs(caption = paste0("correlation ", cor(df$ctg, df$dapi, use = "complete.obs") %>% round(digits = 4))) + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ctg_dapi_well_perline.pdf"), width = 6, height = 6)

df %>%
  mutate(line = substr(id, 1, 7)) %>%
  ggplot(aes(actin, ctg)) + 
  geom_point() + 
  geom_hline(yintercept = 1) + 
  geom_smooth(se = FALSE) +
  theme_cowplot() + 
  facet_wrap(~line) + 
  labs(caption = paste0("correlation ", cor(df$ctg, df$actin, use = "complete.obs") %>% round(digits = 4))) + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ctg_actin_well_perline.pdf"), width = 6, height = 6)

df %>%
  mutate(line = substr(id, 1, 7)) %>%
  ggplot(aes(permeability, ctg)) + 
  geom_point() + 
  geom_hline(yintercept = 1) + 
  geom_smooth(se = FALSE) +
  theme_cowplot() + 
  facet_wrap(~line) + 
  labs(caption = paste0("correlation ", cor(df$ctg, df$permeability, use = "complete.obs") %>% round(digits = 4))) + 
  ggsave(paste0(PATH, "reports/figures/drug_viability/gg_ctg_perm_well_perline.pdf"), width = 6, height = 6)

cor_df <- df %>% 
  mutate(line = substr(id, 1, 7)) %>% 
  nest(-line) %>% 
  mutate(r = purrr::map(data, ~ .x %>% as.data.frame() %>% 
  column_to_rownames("id") %>%
  cor(use = "complete.obs") %>% head(1) %>% as.data.frame())) %>% 
  unnest(r) %>% 
  as_tibble() %>% 
  dplyr::select(-data) %>% 
  as.data.frame()

cor_df %>% write_csv(here::here("reports/tables/source_data_ctg_correlation_heatmap_2i.csv"))

cor_df %>% 
  dplyr::select(-ctg) %>%
  column_to_rownames("line") %>% 
  pheatmap::pheatmap(cluster_cols = FALSE,
                     filename = here::here("reports/figures/drug_viability/ctg_correlation_heatmap.pdf"))

cystic vs solid organoid lines

#UMAP Cystic (Lines 18, 13, 27, 30) vs. Solid (others) treated with DMSO, for Figure 1 / matching expression analysis done for cystic vs. rest

set.seed(123)

cystic_l <- organoid_morphology %>% filter(morphology == "cystic") %>%.$line %>% paste0(., "01")
dense_l <- organoid_morphology %>% filter(morphology == "solid") %>%.$line %>% paste0(., "01")

df <- umap_df %>%
  filter(drug == "DMSO") %>% 
  filter(partition %in% c(1,2)) %>%
  mutate(morphology = case_when(line %in% cystic_l ~ "cystic",
                                line %in% dense_l ~ "solid",
                                 TRUE ~ "other"))

gg_cystic <- umap_df %>% 
  ggplot(aes(v1, v2)) + 
  geom_point_rast(alpha = 1, size = 0.35, color = "#f1f1f1") + 
  
  geom_point_rast(data = df %>%
    filter(morphology != "other") %>%
      sample_frac(1),
  aes(color = morphology),alpha = .1, size = 0.35, shape=16) +
  
  # geom_density_2d(data = df %>% #geom_density_2d_filled
  #   filter(morphology != "other"), # %>% 
  #   #  sample_frac(0.05)
  #   aes(fill = morphology), size = 1.5) +

  scale_color_brewer(type = "qual") +
  #scale_fill_manual(values = c("#0571b0", "#ca0020")) +    
  #scale_color_manual(values = c("#0571b0", "#ca0020")) +   
  #geom_density2d(color = "black") + 
  theme_classic() +
  labs(x = "UMAP 1",
       y = "UMAP 2")+
       #caption = "control treated organoids") + 
  theme_cowplot(font_size = 8) + 
  #theme(legend.position = "nothing")  + 
  coord_fixed()


gg_cystic + ggsave(here::here("reports/figures/mofa/gg_cystic.pdf"), width = 4, height = 4)

Plot Export

plot_grid(plot_grid(gg_size_dist_morph_ridge, gg_size_replicate, labels = c('A', 'B'), label_size = 12, ncol = 2),
          gg_size,
          plot_grid(gg_line, gg_cystic, labels = c('D', 'E'), label_size = 12, ncol = 2),
          labels = c('', 'C', ''), label_size = 12, ncol = 1) +
  ggsave(paste0(PATH, "reports/panels/imaging/panel_size_dist.pdf"), width = 8, height = 16)

plot_grid(plot_grid(ggdrug_size, ggdrug_count, labels = c('A', 'B'), label_size = 12),
          gg_drug,
          gg_size_drug,
          labels = c('', 'C', 'D'), label_size = 12, ncol = 1) +
  ggsave(paste0(PATH, "reports/panels/imaging/panel_size_drug.pdf"), width = 8, height = 12)

gg_size_supp
sessionInfo()
## R version 4.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggridges_0.5.2  scico_1.1.0     princurve_2.1.4 cowplot_1.0.0  
##  [5] ggrastr_1.0.1   here_0.1        readr_1.3.1     tibble_3.0.1   
##  [9] purrr_0.3.4     magrittr_1.5    tidyr_1.1.0     dplyr_1.0.0    
## [13] ggplot2_3.3.1  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.4.6       RColorBrewer_1.1-2 plyr_1.8.6         vipor_0.4.5       
##  [5] pillar_1.4.4       compiler_4.0.0     tools_4.0.0        digest_0.6.25     
##  [9] lattice_0.20-41    nlme_3.1-147       evaluate_0.14      lifecycle_0.2.0   
## [13] gtable_0.3.0       mgcv_1.8-31        pkgconfig_2.0.3    rlang_0.4.6       
## [17] Matrix_1.2-18      yaml_2.2.1         beeswarm_0.4.0     xfun_0.14         
## [21] withr_2.2.0        stringr_1.4.0      knitr_1.28         generics_0.0.2    
## [25] vctrs_0.3.0        hms_0.5.3          rprojroot_1.3-2    grid_4.0.0        
## [29] tidyselect_1.1.0   glue_1.4.1         R6_2.4.1           ggbeeswarm_0.6.0  
## [33] rmarkdown_2.2      pheatmap_1.0.12    farver_2.0.3       splines_4.0.0     
## [37] codetools_0.2-16   backports_1.1.7    scales_1.1.1       ellipsis_0.3.1    
## [41] htmltools_0.4.0    colorspace_1.4-1   labeling_0.3       stringi_1.4.6     
## [45] munsell_0.5.0      crayon_1.3.4       Cairo_1.5-12